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ABSTRACT 



A series of experiments was carried out with two streamlined 
foils at two different speeds in a towing tank. The ambient 
turbulence (nearly isotropic) was generated by means of a bi- 
planar grid. The model-grid separation distance was varied 
systematically in order to subject the trailing vortices to 
varying degrees of turbulence. The data have been expressed 
in terms of a normalized turbulence parameter and the relative 
rise of the vortices. The results have shown that the ultimate 
height to which the vortices rise prior to their demise is 
controlled by the two parameters cited, irrespective of the 
shape and the aspect ratio of the lifting surfaces. A numeri- 
cal analysis has been initiated to analyze the motion of the 
trailing vortices and the propagation of internal waves in 
density stratified medium with the ultimate objective of incor- 
porating the effect of turbulence into the numerical model. 
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INTRODUCTION 



I . 



Vortices and vortex wakes have become a major theme of 
aerodynamics research since the advent of the large aircraft 
and the understanding of their evolution required an examination 
of many fundamental problems in fluid mechanics. Much of the 
progress made during the past two decades was discussed at 
the Symposium on Aircraft V7ake Turbulence and Its Detection 
[Ref. 1] .and at the Aircraft Wake Vortices Conference (Ref. 2]. 
Comprehensive reviews of the entire subject have been given by 
Donaldson and Bilanin [Ref. 3], VJidnall [Ref. 4] and Hallock 
and Eberle [Ref. 5] . 

These studies, as well as numerous others carried out since 
1977, have uncovered a number of complex problems which must 
be resolved in order to achieve a better understanding of the 
important features of trailing vortices in homogeneous and 
stratified media. The principal ones are as follows [Ref. 6]: 

(a) Roll-up process: The velocity and turbulence distribu- 

tions at any station behind the wing depend on the wing 
section, wing-tip shape, Reynolds number, wing incidence, 
and the distance of the station from the wing [Ref. 7] . 
The distributions of the initial velocity and turbulence, 
which influence the roll-up and the decay process, 
cannot be changed independently. For example, a change 
in tip shape changes the core size, as well as the 
velocity and turbulence distributions. High. levels of 
turbulence result in an increased diffusion of vorticity, 
which in turn increases the core size. 

(b) Probe sensitivity of the vortices: Flow visualization 

studies suggest that trailing vortices are extremely 
sensitive to disturbances created by even very small 
probes or bubbles. This forces one to use non- intrusive 
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means of measurements such as an LDV. Even then, 

'vortex wandering' (Ref. 8], which makes the vortices 
appear larger than normal in time-averaged velocity 
measurements (for vortices generated by a wing in a 
wind .tunnel), or the unsteady nature of the flow 
(for vortices generated by a wing in a tow basin) makes 
the mean velocity profiles in the vortices difficult to 
determine . 

(c) Large-scale instabilities; The vortices are seldom 
observed to decay away owing to viscous and turbulent 
dissipation, but are almost always destroyed by either 
mutual induction instability (Crow instability [Ref. 9]) 
and/or vortex breakdown. The Crow instability grows 
exponentially, and results either in a linking of the 
vortex pair into a series of crude vortex rings or in 

a highly disorganized intermingling of the vortices. 
Vortex breakdown, whose mathematical details have not 
yet been adequately treated, rearranges the vortex 
structure and increases the core size, turbulence, and 
energy dissipation. Thus, it is very difficult to 
measure accurately the trajectories of the three- 
dimensional vortices from their creation to their 
ultimate demise. 

(d) Reynolds number; Even the highest Reynolds numbers, 
based on wing chord, reached in wind tunnels or towing 
basins, are an order of magnitude lower than what is 
possible for an aircraft. Thus, the scale effects are 
not easy to assess. 

(e) Ambient conditions such as turbulence and stratification 
play major roles in the evolution of vortices. The 
quantification of these effects requires numerical 
analysis and extremely careful experiments. 

(f) Ground or free-surface effects; The vortex pair may 
move towards a rigid boundary at which the no-slip 
condition must be satisfied or towards a free surface 

at which the zero-shear condition must be satisfied. In 
either case, the vortices come under the influence of 
their images and move accordingly. 

The phenomenon is further complicated by several additional 
facts. When the vortices are propelled towards a rigid sur- 
face, vorticity of opposite sign is generated on the no-slip 
boundary and swept towards the vortex pair. The total vorticity 
diminishes very quickly as vorticity from the two regions 
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diffuses, the wall region serving as a strong sink for the 
vorticity associated with the original vortex (Ref. 10]. The 
development of a boundary layer along the rigid wall may give 
rise to flow separation for sufficiently high Reynolds numbers. 
With or without such a separation, however, the center of the 
vortex pair eventually moves away, or 'rebounds,' from the 
wall (Refs. 10-12]. 

For the case of a zero-stress boundary, the free surface 
still acts as a vorticity sink, but this is relatively weak 
due to the absence of intense oppositely-signed vorticity. 

Thus, in the absence of other impeding phenomena,- one expects 
a mild interaction between the vortices and the free surface 
and a small rebound of the vortex center from the free sur- 
face. However, the ability of the free surface to deform under 
the influence of strain fields leads to a strong interaction 
between the vortices and the free surface. 

It is evident from the foregoing that the motion and the 
life-span of trailing vortices are governed by a number of 
nonlinearly-dependent complex phenomena. A number of experi- 
mental and analytical studies have been carried out at the 
Naval Postgraduate School by Sarpkaya and his students (Refs. 

6, 13-16] in order to investigate the effects of these param- 
eters on the rise and demise of the trailing vortices in 
homogeneous and density stratified media. These studies have 
clearly identified the various demise mechanisms in both media 
and established basic relationships between the rise of vortices 
and the governing parameters in a finite as well as effectively 
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infinite medium [Ref. 17], free from ambient turbulence. The 
present investigation is a continuation of the foregoing studies 
and is confined to the determination of the effect of ambient 
turbulence on the rise and demise of trailing vortices. 
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II. EXPERIMENTAL EQUIPMENT AND PROCEDURES 



A. EQUIPMENT 

The equipment used to generate the trailing vortices has 
been extensively used at this facility over the past three 
years (Refs. 6, 14-17]. Only the salient features, most recent 
modifications and the adaptation for this work are briefly 
described in the following. 

The system consists essentially of a towing basin. The 
auxiliary components of the basin are the plumbing for water 
and for water and brine mixture (to establish a density strati- 
fication) , turbulence management system, top and bottom carriages, 
velocity measuring system, lighting system, and the models [Refs. 
14-17]. 

A mirror, inclined 45 degrees from the horizontal, is attached 
to the top of the tank at the test section to allow viewing of 
the vortex pair and the surface disturbances from the top of 
the basin. 

The drains are provided at the bottom at each end of the 
basin. Two parallel rails are mounted along the bottom of the 
tank. A carriage rides smoothly on these rails and provides 
the test body with a constant velocity through the use of an 
endless cable and a DC motor. The velocity of the model is 
measured and continuously monitored through the use of a mag- 
netic linear displacement transducer. It yields a signal 
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proportional to the velocity of the model with an error less 
than 1 percent. 

The two rails, the carriage and the filling pipes are located 
on or near the bottom of the basin and under a turbulence manage- 
ment system. The system consists of one inch thick polyurethane 
foam sandwiched between two perforated aluminum plates. 

B. MODELS AND GRIDS 

Two rectangular foils (NACA 0012) were used in the present 
study (RPl with B = 6.8 in., b^ = 5.34 in., and c = 3.50 in., 
and RP2 with B = 4.50 in., b^ = 3.53 in., and c = 2.32 in.). 

The interior of these models was hollowed and used as a dye 
reservoir to seed the vortex cores. The models are mounted on 
their bases by means of a thin streamlined aluminum bar with a 
cross section of a NACA 0006 foil and set at the desired angle 
of attack. As noted earlier, all models are pulled by means of 
a DC motor, pulley, and cable system at the desired speed 
(ranging from 0.8 to 4.0 feet/second). Additional details of 
the model construction and mounting are given in Ref. 13. 

Two grids were constructed but only one was used in the 

experiments. The test grid consisted of square plexiglass 

bars (3/8 in. x 3/8 in.). One layer was placed orthogonally 

on top of another layer so as to form a biplanar grid. The 

mesh size of the grid was M = 2.5 in. (see Fig. 1). The 

2 2 

solidity of the mesh, defined as (2Md-d ) /M was 0.2775. The 
grid was attached to the top carriage and towed at desired 
speeds and x/M distances ahead of the model. It was made sure, 
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through repeated experiments, that both the model and the grid 
ran smoothly, with very little or no velocity fluctuations. 

C. GRID-GENERATED TURBULENCE 

A large number of grid- turbulence experiments have been 

conducted in wind tunnels and in water flumes [Refs. 18-23] . 

The measurements have shown that the grid-generated turbulence 

2 2 2 

is only approximately isotropic with U 2 = u^ = 0.83 for a 
uniform duct where u^^ , U 2 , u^ are the longitudinal, lateral 
and vertical velocity components, respectively. In general, 
the mean-square velocity components of the turbulence decrease 
with the downstream distance according to the power law [Ref. 
21 ] : 



where A^ are constants, M is the mesh size of the grid, U is 
the mean velocity of the ambient flow, and the power n ranges 
from 1.25 to 1.5. A careful examination of the data presented 
in Refs. 18-23 has shown (see Fig. 2) that the most suitable 
values of A^^ and n for the grid used in the present investi- 
gation are A^^ = 0.05 and n = 1.3. Thus, one has 




A^(x/M) ^ (i = 1 to 3) 



( 1 ) 




0.05 (x/M) 



-1.3 



( 2 ) 



Consequently, the twice the turbulent kinetic energy per unit 
mass becomes. 



(3) 



2 / 2 ^ 2 ^ 2 , 
q = +U 2 + 0 ^) 



((0.05 +2 x0.83)U^ (x/M)"^*^ 



The dissipation rate e of the grid-generated turbulence per 
unit mass is defined as 



e = -d/2)dq^/dt 



(4) 



and decreases with increasing x/M. Using Eq. (1) and the 
definitions of e and q, one has 



= (n/2) (x/M)"^ (5) 

q U 

Friehe and Schwarz (Ref. 20] were able to collapse all of their 
grid turbulence data onto a single curve through the use of 
Eq. (4). For n = 1.3, one has 



= 0.65(x/M)~^ (6) 

q U 



and 



^ = 0.0865(x/M)"^'^ (7) 

U 



The turbulence parameter e* , defined by Crow (Ref. 24] as. 



£ * 



(eb^) 




( 8 ) 
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may be combined with Eq. (7) to yield 



1 0.4423 



17 ^ (X/M) 



-0.767 



(9) 



(V^/U) (M/b^) 



in which M/b^ represents the ratio of the mesh size to the 
initial vortex separation. The values of V^/U and M/b^ are 
known for a given grid, model, and model velocity. It must 
also be noted that the exact value of (0.05 in Eq. (2)) is 
not very critical. It is easy to show that a change even as 
large as 20 percent in (say from 0.05 to 0.06) causes only 
a 6 percent change in e*. 

The integral scale of the grid-generated-turbulence 

field has been measured by a number of investigators (see e.g.. 
Refs. 18-23). The results have shown that increases with 

x/M. Huot et al. [Ref. 23] found experimentally that 
m.ay be represented by 



in which x^/M represents the distance to a fictitious origin. 
For the grid size used in the present experiments, x^/M is 
about 2.5. However, it is not necessary to retain x^/M in 
Eq. (10) simply because in the present experiments x/M ranged 
from a minimum value of 50 to about 1000. Thus, one has 



Lii/M = 0.14 (x/M - x^/M) 



0.4 



( 10 ) 



Lii/bo = 0.14 (M/b^) (x/M) 



0.4 



( 11 ) 
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Equation (11) shows that is about 1.57 for the smaller 

model (RP-2 with M/b^ = 0.707) and about 1.0 for the larger 
model (RP-1 with M/b^ = 0.47). Thus, it is seen that the inte- 
gral scale of the turbulence generated by the grid used in the 
present experiments varies from one half to 1.5 times the 
initial vortex separation. In other words, it is within the 
distances cited that the integral scale of the background in- 
creases and becomes comparable to the size of the individual 
vortices . 

D. TEST PROCEDURES 

A model is placed in the basin and the basin is filled 
gradually with fresh water to the desired level. In the cases 
where dye was used to seed the vortex pair, the models were 
filled with dye prior to filling the tank. After removal of 
trapped air and after a sufficient period of waiting for the 
escape of dissolved air in the water and the elimination of 
any internal currents in the basin, the grid was set in motion 
at the desired speed. When the grid has moved a prescribed 
distance of x/M, the model was set in motion. As noted earlier, 
the grid and the model moved at identical speeds, i.e., the 
horizontal distance between them was kept constant at a given 
value of x/M. Most of the experiments were repeated at least 
three times. 

The motion of the trailing vortices are recorded on high- 
speed Polaroid film at the test section (one of the plexiglass 
panels near the middle of the basin) . Each picture included 
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two clocks accurate to 0.1 seconds, the vertical and horizontal 
scales on the plexiglass window and, of course, the side view 
of the trailing vortices as they rose from the model after 
formation. The time interval between successive pictures is 
determined from the two clocks. The first picture always cap- 
tured the instant the grid arrived at the test section. The 
second picture captured the instant the model passed through 
the test section. The subsequent pictures (taken at about 0.75 
second intervals) captured the rise and demise of the vortices. 
The vertical rise of the vortices is determined from the ver- 
tical scale. Attention has been paid to the fact that the 
vortices are farther away from the scale on the window and 
that the scale placed vertically in the middle of the test 
section does not exactly correspond to the scale marked on the 
window due to refraction and parallax. The necessary correc- 
tion was made by photographing a scale placed in water in the 
middle of the test section together with the scale marked on 
the window. This resulted in a simple conversion table which 
enabled the determination of the actual position of the vortex 
from the scale reading on the photograph. 

The results are normalized and plotted in various forms 
and compared with those obtained in the previous runs. Each 
experiment was repeated three times for most of the data pre- 
sented herein. All trailing vortices were recorded on film 
until the time they have completely dissipated either due to 
aging (diffusion of vorticity due to viscosity, turbulence, 
entrainment, and detrainment) , or due to Crow instability 
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(sinusoidal instability and linking of vortices), and/or due 
to vortex breakdown (core bursting) . It was thus possible to 
determine the life span of vortex cores from their formation 
to their final demise. 
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III. ANALYSIS 



A. DIMENSIONAL ANALYSIS 

The dependent parameter of major importance is the instan- 
taneous position of the vortex pair (x,y). It may be expressed 
as a function of the following parameters [Ref. 17]: 

X = f {t,U,d^,p^,dp/dy,v,B,Ar,a,g,r^,e,L^j^) (12) 

and 

y = f (t,U,d^,p^,dp/dy,v,B,AR,a,g,r^,e,Lj^^) (13) 

in which t represents the time, U the velocity of the model, 
d^ the initial depth of the vortex pair, p^ the reference 
density of water, dp/dy the linear density gradient, v the 
kinematic viscosity of water, B the base width of the model, 

AR the aspect ratio of the model, a the angle of attack of the 
model, g the gravitational acceleration, a the surface tension 
of water, r^ an effective core radius, characterizing the 
effect of the wing-tip shape in addition to other wing param- 
eters, e the rate of decay of the turbulent energy per unit 
mass, and L^^j^ the integral scale of the turbulent field. 

The height and width of the test- section and the height of 
the model within it were not included in the foregoing because 



23 



a detailed analysis, based on ideal vortices, has shown that 
the velocities induced by the bottom were negligible. Effects 
of surface tension on the instantaneous position of the vortex 
pair is deemed negligible. It is for this reason that o is 
not included in Eqs . (12) and (13). 

A dimensional analysis of Eqs. (12) and (13) yields 

x/B = f (Ut/B,d^/B,NB/U,U^/gB,UB/v,AR,a,r^/B,e*,Lj^^/b^) 

(14) 

y/B = f (Ut/B,d^/B,NB/U,U^/gB,UB/v,AR,a,re/B,e*,Lj^j^/b^) 
in which 



N = (-g/p^ dp/dy)^/^ (15) 

and is known as the Brunt-Vaisala frequency. The other varia- 
bles have been defined previously. 

Equation (14) may be recast in terms of the initial separa- 
tion b of the vortex pair, and the initial mutual induction 
velocity of the vortices by noting that B/b^ and V^/U are 
uniquely determined by AR and a for a given wing shape. Thus, 
one has 



x/b = f (V t/b ,d /b ,Nb /V ,V /gb ,V b /v,r /b ,e*,L, ,/b ) 

' o o ' o' o' o o' o ' o' ^ o ' o o' e o 11' o 

( 16 ) 



y/b = f (V t/b ,d /b ,Nb /V ,V"^/gb ,V b /v,r /b ,e*,L, ,/b ) 
o o ' o o o 0 0 0 ^ 0 0 0 e o 11 o 
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Experiments have shown that [Ref. 6] the ratio of the 
initial vortex-core spacing to wing span, b^/B, is 0.70 for the 
sharp-edged Delta wing, 0.73 for the round-edged Delta wing, 
and nearly equal to tt/ 4 for the rectangular wings (NACA 0012) . 
Whereas all the parameters in Eq. (16) may be changed indepen- 
dently, is taken as nature provides it. The primary 

reason for this is that a century of theoretical and experi- 
mental aerodynamics research has been incapable of describing 
the details of the structure of the tip vortex to be used as 
initial conditions in a viscous solution. It is surprising, 
but true, that until recently the importance of the wing-tip 
shape and its influence upon both the initial tangential velocity 
profile and the initial turbulence in the vortex has not been 
fully appreciated. Here the said influence has been charac- 
terized in terms of an effective core radius with full awareness 
of its shortcomings. 

B. NUMERICAL ANALYSIS OF THE VORTEX MOTIONS 

The generation of the internal waves and the rise and demise 
of the vortices in a stratified medium may be analyzed through 
the use of the equations of motion for an incompressible fluid 
for both laminar and turbulent motions provided that a suitable 
turbulence closure model is adopted and the usual Boussinesq 
approximation (gravitational acceleration is much larger than 
the fluid accelerations) is made. Fgr the type of motions con- 
sidered herein the Boussinesq approximation is quite valid and 
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has been used in the investigation of all types of internal 
waves in stratified fluids. 

For a two-dinensional flow, with y vertical and x horizon- 
tal, the equations of motion are 
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and the equation of continuity 
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Defining vorticity as usual by 
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and eliminating pressure between Eq. (17) and Eq. (18), one 
has 
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in which p ' 



is the fluctuating part of the density p given by 



p = p^ + p (y) + p ' (x,y , t) 



( 22 ) 
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where is the reference density, and p(y) the initial den- 
sity at elevation y at t = 0 . 

The diffusion of density is given by 



Combining Eqs. (19) , (22) and (23) , and simplifying, one has 



It is convenient to cast the foregoing equations into non- 
dimensional forms, scaling each variable by a quantity charac- 
teristic of its expected magnitude. Two possible time scales 
exist. There is the dynamic time scale, which is the time a 
characteristic length would be traversed by a fluid particle 
traveling at the characteristic velocity, and there is the 
buoyant time scale based on the natural buoyancy frequency of 
the stratified flow, i.e., the Brunt-Vaisala frequency N, 
defined by Eq. (15) . Each of the scales gives a slightly 
different form of the normalized governing equations. 

Dynamic Scale: 

Introducing and as the characteristic velocity and 
length, one has 
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Normalizing Eqs . (21) and (24), one has 
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in which 
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Equations (25) and (26) are valid when >> 1, i.e., when 
the buoyancy has very little influence on the nonlinear dynamics 
of the motion. In this case, the density perturbation acts as 
a passive scalar advected by the velocity field. 

Buoyant Scale : 

Introducing the following dimensionless parameters: 
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and normalizing Eqs. (21) and (24) , one has 





( 29 ) 



and 
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F n V 
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(30) 



2 

in which Re and n are as defined in Eq, (27) 



Equations (29) and (30) are valid when F^ << 1 and buoyancy 
dominates the flow. When approaches zero, the equations 
that result from Eqs. (29) and (30) describe the propagation 
of linear internal waves. The F^ << 1 regime is of interest 
in the present investigation because for submerged bodies of 
naval interest F^ (the Froude number) is about 0.01. The 
analysis will consider the full nonlinear equations given by 
Eqs. (29) and (30) rather than their linearized form (i.e., 
the case of = 0). However, the flow will be assumed to be 
inviscid, i.e., the viscous diffusion will be ignored to a 
first order approximation. Subsequent analysis will incorpor- 
ate the effects of viscous and turbulent diffusion into the 
numberical calculations. The velocities u and v are given by 
the Biot-Savart law as 
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and 



v(x,y) = / dx’dy' 

2irr'^ 



(32) 



in which x,y is the point at which the velocity is calculated 
and x',y* is an arbitrary point at which the vorticity is 
C (x ' ,y ' ) . 

The radial distance r is given by; 



r^ = (X-X')^ + (y-y')^ (33) 

Normalizing Eqs . (31-33) through the use of the charac- 
teristic dimensions given in Eq . (28), one has: 
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It is appropriate to take v^ , the initial mutual induction 
velocity of the vortex pair, as the characteristic velocity 
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and the , the initial separation of the vortices, as the 
characteristic length. Thus, one has: 
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In dimensionless forms, one has: 
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in which N^b^/v^ is called the stratification parameter SP. 

The value of SP normally varies from zero to about 100. The 
larger the SP, the stronger the stratification. For example, 
for a model running at a speed of 5 ft/sec at an angle of 
attack of 10 degrees, SP smaller than about 1.0 represents weak 
stratification, SP between about 1 and 5 represents medium 
stratification, and SP larger than 5 corresponds to strong 
stratification . 

The numerical solution of Eqs. (29) and (30) through the 
use of Eqs. (34) -(36) and the appropriate boundary conditions 
will be discussed later. 
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IV. DISCUSSION OF RESULTS 



A. EFFECTS OF AMBIENT TURBULENCE ON VORTEX MIGRATION 

Figures 3-8 show the trajectories of the vortex cores for 
the model RPl, moving at U = 0.9 ft/sec, for various values of 
X/M. In each figure H* = H/b^ and T* = V^t/b^. The solid line 
represents H* = T* which is the theoretical result [Ref. 6] 
for the migration of a pair of vortices in a smooth, inviscid, 
homogeneous medium. Figures 9-14 show similar results for the 
same model at a speed of U = 1.6 ft/sec. Figures 15-22 show 
the vortex trajectories for the model RP2 , moving at U = 0.9 
ft/sec. Finally, Figures 23-31 show similar data for RP2 at 
a model speed of U = 1.6 ft/sec. 

Figures 3-31 show that there is some scatter in the data. 
The maximum migration of vortices, at which time they are 
completely destroyed by instability events, decreases with 
increasing ambient turbulence. The scatter is partly due to 
the random nature of the ambient turbulence, and partly due to 
the fact that the vortex pair does not deform or demise in a 
similar manner either in one test or from one test to another. 
The three-dimensional nature of the ambient turbulence, the 
rotation of the pair along a horizontal axis (under certain 
circumstances), sinusoidal instability, and the occurrence of 
vortex breakdown (in only one of the vortex pair) change the 
relative positions of the vortices at the test section and 
give rise to scatter in the data. It is remarkable that the 
scatter is not any larger. 
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Figures 32-35 are composite plots of Figures 3-8, 9-14, 
15-22 and 23-31, respectively. This same data are presented 
in tabular form in Appendices B-D in terms of T* , H* and e*. 

For the grid used in the tests, the normalized integral 
scale of the turbulence, is given in Table I for the 

X/M values encountered in the tests. 

TABLE I 

NORMALIZED INTEGRAL SCALE OF TURBULENCE 
(Lj^j^/b^) FOR VARIOUS X/M VALUES 

x/M ^ll/^o ^ll'^^O 



00 



1000 


1.04 


1.57 


750 


0.93 


1.40 


500 


0.79 


1.19 


250 


0.60 


0.90 


180 


0.52 


0.79 


80 


0.38 


0.57 


50 


0.31 


0.47 


25 


0.24 


0.36 



The integral scale of the turbulence field for RPl is comparable 
to the separation distance of the vortices, for x/M larger than 
about 500. The same is true for RP2 for x/M larger than about 
180. Figures 32-34 show that the path of the vortices in the 
turbulent field does not significantly deviate from the path of 
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the vortices in the non-turbulent medium until the vortex pair 
rises at least one initial separation distance. Subsequently, 
the vortices slow down considerably before the onset of insta- 
bility and eventual dissipation. It is during this decelera- 
tion period that the integral scale becomes important in 
bringing about the instability which breaks up the vortices. 

In other words, strong turbulence with small integral scale 
(small x/M) leads to the eventual deposition of the circulation 
in larger and larger regions of the vortex without the more 
dramatic and destructive effects of larger scale motions. On 
the other hand, weaker turbulence with larger integral scale 
(for large x/M) leads to strong instabilities, breaks up the 
vortices and leaves them in a state to be devoured by small 
intensity turbulence. Subsequently, the vorticity is diffused 
outward and finally destroyed completely. It is clear from the 
foregoing that the lifespan of the vortices depends significantly 
both on the scale and the intensity of the ambient turbulence. 

The quantification of the separate roles of these two param- 
eters needs additional experiments with other grids. 

Figures 32-34 also show that the smaller x/M, the smaller 
is the ultimate rise of vortices. Leaving aside for the time 
being the rise history of the vortices and the separate effects 
of the integral scale and turbulence intensity, the attention 
will be focused now on the maximum height attained by the vor- 
tices. The said height is of extreme importance in practical 
situations. Figure 36 shows the maximum height data obtained 
with both models as a function of e*. Also shown in this figure 
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is the flight data obtained by Tombach (Ref. 25] . There is 
reasonable agreement between the laboratory and flight data 
in spite of the considerable scatter due to the difficulty of 
determining the maximum height attained by the vortices. As 
noted earlier, the uncertainties stem partly from the random- 
ness in the occurrence of the instabilities leading to the 
destruction of the vortices and partly from the experimental 
difficulties and uncertainties (diffusion of dye, subjectivity 
in deciding the exact position of the vortex core, etc.). Evi- 
dently, the turbulence parameter is the primary governing 
parameter in the determination of the maximum height. The 
role of the integral scale is not clearly discernible and is 
partly obscured by the scatter in the data. 

Figure 37 shows the dimensionless time t* {= ^Q^/b^) as a 
function of the turbulence parameter. The data points shown 
on the T* axis correspond to the no-turbulence case. It is 
assumed that the effect of turbulence is negligible for e* 
smaller than about 0.1 in order to show the entire data in the 
figure. Clearly, the effect of turbulence is to reduce both 
the lifespan and the maximum height attained by the vortices. 
The data have also shown that the vortices in the non-turbulent 
medium break up primarily due to linking [Ref. 6] and mostly 
due to vortex breakdown in the turbulent medium. Thus, it is 
concluded from Figures 36-37 and the numerous observations that 
turbulence enhances the instabilities leading to the vortex 
breakdown and, at the same time, diffuses vorticity rapidly. 
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eventually leading to the total destruction of the vortices. 
This mechanism is far more powerful than the sinusoidal insta- 
bility leading to the linking of the vortices. It remains to 
be determined as to what the individual effects of the turbu- 
lence intensity and scale are in bringing about the destruction 
of the vortices due to the vortex breakdown. 

B. EFFECTS OF STRATIFICATION ON VORTEX MIGRATION 

Equations (29) and (30) were integrated through the use of 
an efficient upwind differencing scheme [Ref. 26] for the 
inviscid flow case, i.e.. Re = ®. A copy of the complete com- 
puter code is given in Appendix E. 

The initial vorticity distribution was assumed to be 
Gaussian , i.e.,: 
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^ ^ m' om ^ 



(39) 



The calculations had to be carried out in a finite mesh, 
as in all other numerical analysis. The approximate boundary 
conditions are shown on the boundaries of the quadrant in 
which the calculations have been performed, as well as a 
schematic of the computational nodal points (see Figs. 38a and 
38b) . It should be noted that the boundary conditions 
u(0,y) = v(x,0) = 0 on the axes are automatically satisfied 
by evaluating the Biot-Savart equations in all four quadrants. 
The fields of velocity, vorticity, density (p(y)), and the 
fluctuating components of the density (p') are calculated at 
each time step and plotted at regular intervals. 
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Among the several calculations carried out only one uti- 
lizing the initial parameters corresponding to those of an 
experiment will be reported here. The computer code is quite 
general and can be applied to any stratified and unstratified 
flow situation. As the velocity field is computed from ana- 
lytic expressions (i.e., the Biot-Savart equations) vice an 
algorithm involving an iterative convergence scheme, stability 
of the program is enhanced at a slight expense of execution 
time . 

The plots of the velocity field, constant vorticity, den- 
sity and density fluctuation are shown in Figures 39a-39d 
through Figures 46a-46d for various values of the normalized 
time T*. In the calculations, the initial depth of the 

vortices from the free surface was taken d^/b = 8, the initial 

o o 

value of r /b = 0.09, the stratification parameter Nb /V = 0.75 
o o o o 

and F^ = 0.018. These values correspond to those of experi- 
ments with which the numerical calculations are to be compared 
to. 

Figures 39-46 show that at small times the vortices rise 
vertically and the circulation in the flow field is primarily 
due to the initial vortices. As time increases (see e.g.. 

Figure 43a), the regions of circulation with countersigned 
vorticity develop in the upper right and lower left, regions 
of the vortex. It is easy to show that this counter vorticity 
not only reduces the rise velocity of the original vortex 
pair, but pushes the pair against each other. Consequently, 
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the vorticity is lost in the overlapping regions of the vortex 
pair, and the rise velocity is further reduced. As time in- 
creases, the countersigned vorticity begins to dominate the 
flow (see e.g., Figure 44a) and the vortex migration stops. 

With further increases in time the vortex begins to migrate 
downward . 

The density contours reveal the same phenomenon in a dif- 
ferent context. As the vortices rise, fluid of greater den- 
sity is pushed upwards (see e.g., Figure 43c) into regions of 
lesser density. Since such a migration cannot go on indefinitely, 
the vortices rise to a maximum height and then begin to sin)c 
downwards. The calculations do not take into account the sinu- 
soidal instability and the vortex breakdown. Consequently, 
the vortices in the numerical calculations continue to exist 
unimpeded by the instability mechanisms. In reality, the 
vortices begin to break up as they near the end of their maxi- 
mum migration and eventually disappear. 

The experimental and calculated values are compared in 
Figure 47. The correspondence between the measured and calcu- 
lated values is surprisingly good up to the time of maximum 
rise. This is partly because of the experimental fact that the 
rise of vortices in a highly stratified medium (Nb^/V^ = 0.75, 

F = 0.018) is not strongly dominated by sinusoidal instability 
or vortex breakdown. The migration of vortices is inhibited 
primarily due to the reduction of vorticity of the initial 
vortices and the creation of countersigned vorticity. 
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Additional calculations are underway to examine the 
effects of viscous diffusion and ambient turbulence for 
various degrees of initial stratification. The results 
presented herein and the initial comparisons are extremely 
encouraging and are expected to lead to a better understanding 
of this extremely complex and challenging phenomenon. 
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V. CONCLUSIONS 



The investigation reported herein warranted the following 
conclusions : 

1. The effect of ambient turbulence in homogeneous medium 
is to inhibit the rise of the vortices; 

2. The larger the intensity of the ambient turbulence, 
the smaller is the ultimate rise of vortices; 

3. The phenomenon seems to be governed primarily by the 
turbulence parameter e*. The effect of the integral 
scale of the turbulence field seems to be secondary, 
at least for integral scales on the order of initial 
vortex separation; 

4 . The demise mechanism appears to be the precipitation of 
the instabilities (primarily vortex breakdown) by turbu- 
lence, the breaking up of the vortex filaments, and the 
eventual destruction of the remaining vorticity by turbu- 
lent diffusion; 

5. Numerical calculations regarding the migration of vortices 
in a non-turbulent stratified medium have shown that the 
effect of stratification is to reduce the rate of rise 

of vortices and the maximum height attained by them; 

6. The reduction in maximum height is brought about partly 
by the development of countersign vorticity and partly 
by the destruction of vorticity in the overlapping 
regions of the vortex pair; 

7. The predictions of the numerical model are in good 
agreement with the experimental data. The theoretical 
model could be made more realistic through the inclusion 
of at least the effects of turbulence and possibly those 
of the demise mechanisms. 
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APPENDIX A 



FIGURES 




Figure 1, 



Diagram of Biplanar Grid Construction 
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Figure 2. Comparison of Turbulence Models 
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Figure 3. Vortex Rise for RPl (U = 0.9 ft/sec, X/M 
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Figure 4. Vortex Rise for RPl (U = 0.9 ft/sec, X/M - 750) 
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Figure 5. Vortex Rise for RPl (U = 0.9 ft/sec, X/M = 500) 
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Figure 6. Vortex Rise for RPl (U = 0.9 ft/sec, X/M = 250) 



(0 




in 






50 



Figure 7. Vortex Rise for. RPl (U = 0.9 ft/sec, X/M = 180) 
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Figure 8. Vortex Rise for RPl (U = 0.9 ft/sec, X/M 
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Figure 9. Vortex Rise for RPl (U = 1.6 ft/sec, X/M 
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Figure 10. Vortex Rise for. RPl (U = 1.6 ft/sec, X/M = 750) 
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Figure 11. Vortex Rise for.RPl (U = 1.6 ft/sec, X/M = 500) 
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Figure 12. Vortex Rise for RPl (U = 1.6 ft/sec, X/M = 250) 



U) 




m 



o 



56 



Figure 13. Vortex Rise for RPl (U = 1.6 ft/sec, X/M - 180) 
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Figure 14. Vortex Rise for RPl (U = 1.6 ft/sec, X/M 
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Figure 15. Vortex Rise for RP2 (U = 0.9 ft/sec, X/M 
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Figure 16. Vortex Rise for. RP2 (U = 0.9 ft/sec# X/M - 750) 
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Figure 17. Vortex Rise for RP2 (U = 0.9 ft/sec, X/M = 500) 
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Figure 18. Vortex Rise for RP2 (U = 0.9 ft/sec, X/M - 250) 
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Figure 19. Vortex Rise for.RP2 (U = 0.9 ft/sec, X/M = 180) 
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Figure 20. Vortex Rise for. RP2 (U = 0.9 ft/sec, X/M 
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Figure 21. Vortex Rise for RP2 (U = 0.9 ft/sec, X/M 
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Figure 22. Vortex Rise for RP2 (U = 0.9 ft/sec, X/M 
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Figure 23. Vortex Rise for RP2 (U = 1.6 ft/sec, X/M 
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Figure 24. Vortex Rise for RP2 (U = 1.6 ft/sec, X/M - 1000) 
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Figure 25. Vortex Rise for RP2 (U = 1.6 ft/sec, X/M = 750) 
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Figure 26. Vortex Rise for RP2 (U = 1.6 ft/sec, X/M = 500) 
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Figure 27. Vortex Rise for RP2 (U = 1.6 ft/sec, X/M = 250) 
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Figure 28. Vortex Rise for RP2 (U = 1.6 ft/sec, X/M = 180) 
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Figure 29. Vortex Rise for. RP2 (U = 1.6 ft/sec, X/M 
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Figure 30. Vortex Rise for RP2 (U = 1.6 ft/sec, X/M 
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Figure 31. Vortex Rise for RP2 (U = 1.6 ft/sec, X/M 
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Figure 32a. Comparison of Vortex Rise for RPl for Various 
X/M, U = 0.9 ft/sec, Data Points 
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Figure 32b. Comparison of Vortex Rise for RPl for Various 
X/M, U = 0.9 ft/sec, Curves 
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Figure 33a. Comparison of Vortex Rise for RPl for Various 
X/M, U = 1.6 ft/sec, Data Points 
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Figure 33b. Comparison of Vortex Rise for RPl for Various 
X/M, U = 1.6 ft/sec, Curves 
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Figure 34a. Comparison of Vortex Rise for RP2 for Various 
X/M, U = 0.9 ft/sec. Data Points 
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Figure 34b. Comparison of Vortex Rise for RP2 for Various 
X/M, U = 0.9 ft/sec, Curves 
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Figure 35a. Comparison of Vortex Rise for RP2 for Various 
X/M, U = 1.6 ft/sec, Data Points 
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Figure 35b. Comparison of Vortex Rise for RP2 for Various 
X/M, U = 1.6 ft/sec, Curves 
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Figure 36. Comparison of K* vs. e* for RPl, RP2 at 
ft/sec and 1.6 ft/sec with Tombach Data 
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Figure 37. Comparison of T* ,vs. e* for RPl, RP2 at 
ft/sec and 1.6 ft/sec with Tombach Data 
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Figure 38a. Numerical Coordinate Axis and Boundary 
Conditions 
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Figure 38b. Numerical Nodal Point Arrangement 
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Figure 39a. Velocity Field for T* = 0.0365 
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Figure 39b. Vorticity Field for T* = 0.0365 



88 



Y/BO 



O 



‘ 0 . 0002 • 



• 0 . 0002 - 



CM 



■0.0004- 



■0.0004- 



CO - 



-0.0006- 



•0.0006' 



^ - 



• 0.0008 • 



■ 0 . 0008 • 



LO - 



■ 0 . 0010 - 



■ 0 . 0010 - 



CO - 






GO - 



CD - 



■ 0 . 0012 ' 



•0.0014- 



-0.0016- 



■ 0 . 0012 - 



-0.0014- 



-0.0016- 



-T 1 ^ 1 ^ ^ ^ \ ^ 1 ^ 1 ^ ^ 1 1 1 1 P- 

-5 -4 -3 -2 -1 0 1 2 3 - 4 5 

X/BO 



Figure 39c. Constant Density Contours for T* = 0.0365 
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Figure 39d. Density Perturbation Contours for T* = 
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Figure 40a. Velocity Field for T* = 0.728 
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Figure 40b. Vorticity Field for T* = 0.728 
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Figure 40c. 



Constant Density Contours for T* 



0.728 



93 



Y/BO 



#• 




Figure 40d. Density Perturbation Contours for T* = 0.728 
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Figure 41a. Velocity Field for T* = 1.46 
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Figure 41b. Vorticity Field for T* = 1.46 
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Figure 41c. Constant Density Contours for T* = 1.46 
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Figure 41d. Density Perturbation Contours for T* = 1.46 
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Figure 42a. Velocity Field for T* = 2.18 
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Figure 42b. Vorticity Field for T* = 2.18 
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Figure 42c. Constant Density Contours for T* = 2.18 
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Figure 42d. Density Perturbation Contours for T* = 2.18 
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Figure 43a. Velocity Field for T* = 2.9 
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Figure 43b. Vorticity Field for T* = 2.9 
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Figure 43c. Constant Density Contours for T* = 2.9 
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Figure 43d. Density Perturbation Contours for T'* = 2.9 
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Figure 44a. Velocity Field for T* = 3.64 
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Figure 44b. Vorticity Field for T* = 3.64 
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Figure 44c. Constant Density Contours for T* = 3.64 
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Figure 44d. Density Perturbation Contours for T* = 3.64 
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Figure 45a. Velocity Field for T* = 4.4 
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Figure 45b. Vorticity Field for T* = 4.4 
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Figure 45c, Constant Density Contours for T* = 4.4 
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Figure 45d. Density Perturbation Contours for T* = 4.4 
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Figure 46a. Velocity Field for T* = 5.1 
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Figure 46b. Vorticity Field for T* = 5.1 
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Figure 46c. Constant Density Contours for T* = 5.1 
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Figure 46d. Density Perturbation Contours for T* = 5.1 
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Fiaure 47. Comparison of Numerical and Experimental Results 
for SP = 0.75, F = 0.018 



APPENDIX B 



RPl TABULATED DATA 



Model : RPl 

bQ = 5.34 in 
Us 1.6 ft/s 
ALPHA s 10 Deg 
Vo s 0.030 ft/s 



EPS 1 LON* 


T* ’ 


EPSILON* 


H* 


0.01 


8.2 


0.01 


4.7 


0.01 


5.4 


0.01 


3. 15 


0. 197 


3. 17 


0. 197 


2.5 


0. 198 


3.2 


0. 190 


2.2 


0. 199 


3.0 


0. 199 


2.5 


0.273 


1 . 74 


0.273 


1 . 25 


0.273 


2.2 


0.273 


1 .6 


0.273 


2.4 


0.273 


1 . 04 


0.451 


1 . 74 


0.451 


1 . 43 


0.459 


1 . 72 


0.459 


1.13 


0.595 


0.893 


0.595 


0.89 


0. 595 


1 .07 


0.595 


1 .07 


0. 58 


1 . 27 


0.50 


1 . 25 


1.12 


1 .9 


1.12 


1 .6 


1.03 


1 .0 


1 .03 


1 . 55 


1.00 


1 . 4 


1 .00 


1.17 
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Model; RP1 

bo = 5.34 in 
U = 0.9 ft/s 
ALPHA s 10 Deg 
Vq s 0.018 ft/s 



EPSILON* 


T* 


0.01 


5.5 


0.01 


4. 1 


0. 1 78 


3.3 


0.24 


2.55 


0.24 


2.45 


0.413 


1 . 27 


0.413 


1 . 8 


0. 53 


1 . 4 


0.53 


0.8 


0.97 


1 . 73 



EPSILON* 


H* 


0.01 


3.74 


0.01 


3.15 


0. 178 


2.02 


0.24 


1 . 73 


0.24 


1 . 84 


0.413 


1 . 25 


0.413 


1 .4 


0. 53 


1 .02 


0.53 


0.67 


0.97 


1 . 33 
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APPENDIX C 



RP 2 TABULATED DATA 



Mode 1 : RP2 

— 3 • S 3 in 
U = 0.9 ft/3 
ALPHA = 10 Deg 
V <5 = 0.021 ft/s 



EPSILON* T* 



0. 


01 


4 . 


3 


0. 


01 


3. 


44 


0. 


01 


5. 


4 


0. 


01 


4 . 


1 


0. 


01 


3. 


0 


0. 


1 34 


2. 


0 


0. 


1 34 


4 . 


5 


0. 


1 8 


2. 


1 


0. 


1 6 


1 , 


, 0 


0. 


314 


2. 


, 6 


0. 


319 


1 , 


, 7 


0. 


314 


2, 


2 


0. 


319 


2, 


, 04 


0. 


, 392 


1 , 


. 22 


0. 


39 


1 , 


, 73 


0. 


, 404 


2, 


, 22 


0. 


, 422 


1 , 


, 25 


0, 


, 745 


1 


.23 


0. 


, 75 


1 


. 6 


0, 


, 745 


2, 


. 0 


0, 


, 745 


1 


.65 


1 , 


. 03 


1 


. 5 


1 , 


. 06 


0, 


. 03 


0. 


.97 


1 


. 9 


1 


.0 


1 


. 1 1 


2, 


.02 


0 


. 9 



EPSILON* 


iiH 


0.01 


3.16 


0.01 


3.4 


0.01 


3. 5 


0.01 


3. 1 


0.01 


2.5 


0.134 


2.37 


0. 1 34 


3.05 


0.10 


1 . 5 


0. 10 


1 . 2 


0.314 


1 . 9 


0.319 


1 . 7 


0.314 


1 . 55 


0.319 


1 . 90 


O. 392 


0.97 


0. 39 


1.15 


0. 404 


1 . 6 


0.422 


1 . 2 


0. 745 


1 . 1 


0. 75 


1 . 4 


0. 745 


1 . 9 


0.745 


1 . 7 


1 .03 


1 . 55 


1 . 06 


0.9 


0.97 


1 . 0 


1 .0 


1 . 1 


2.02 


0.9 
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Model: RP2 

&Q = 3.53 in 
U = 1.6 ft/s 
ALPHA = 12 Deg 
Vq = 0.04 ft/5 



EPS 1 LON" 


11 


EPSILON" 


H" 


0.01 


4 . 5 


0.01 


3.91 


0.01 


8 . 4 


0.01 


5.29 


0.01 


5.6 


0.01 


4. 74 


0.099 


3.8 


0.099 


3.48 


0.123 


3 . 35 


0. 1 23 


2. 59 


0. 123 


3.89 


0. 123 


3. 3 


0. 168 


4 . 2 


0. 168 


2.96 


0. 1 69 


5. 1 


0. 169 


3.48 


0.279 


2.97 


0. 279 


1 . 9 


0.284 


3.85 


0.284 


1 . 87 


0. 36 


2 . 86 


0. 36 


1 . 9 


0.366 


2.2 


0. 366 


1 . 87 


0.692 


2.12 


0.692 


1 . 9 


0.682 


1 . 66 


0.682 


1 . 72 


0.692 


1 . 6 


0.692 


1 . 53 


0.945 


1 . 48 


0.945 


1 . 35 


0.989 


1 . 77 


0.989 


1 .62 


0.989 


1 . 35 


0.989 


1 . 6 


1 . 66 


1 . 55 


1 .66 


1.51 



123 



APPENDIX D 



TOMBACH DATA 



EXPERIMENTAL DATA COMPILED FROM DATA REPORTED BY TOMBACH 

( SEE [REF. E4] ) 



EPS 1 LON" 


11 


EPSILON" 


H" 


0.040 


4. 1 


0.040 


4.2 


0.048 


3.5 


0.052 


3.5 


0.055 


3.4 


0.060 


3.4 


0.055 


2.9 


0.059 


2.8 


0.055 


2.4 


0.057 


2. 3 


0.160 


3.8 


0. 165 


3.7 


0.220 


2.2 


0.240 


2.15 


0.220 


1 .6 


0.350 


2.0 


0.220 


1 . 4 


0.420 


1 . 95 


0.170 


1 . 4 


0. 320 


1 . 75 


0.330 


2.0 


0. 390 


1 . 7 


0.400 


2.0 


0.265 


1 . 6 


0. 320 


1 . 8 


0. 350 


1 . 5 


0. 370 


1 . 75 


0. 165 


1 . 3 


0.330 


1 . 55 


0. 265 


1 . 3 


0. 330 


1.2 


0.350 


1 . 2 
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APPENDIX E 



NUMERICAL ANALYSIS PROGRAM LISTING 



PROGRAM V0R1SPL 



5 MARCH 86 



PURPOSE: 

THIS PROGRAM NUMERICALLY SIMULATES THE RISE OF A TRAILING 
VORTEX PAIR SHED OFF A LIFTING SURFACE SUBMERGED IN A 
DENSITY STRATIFIED MEDIUM. 



TO RUN THIS PROGRAM: 

1) GO TO "INPLfT PARAMETERS" SECTION OF CODE AND CHANGE ANY PROBLEM 
PARAMETERS ( SP, FV, DT, ZTMAX, MP, NP, ISAVE) DESIRED. 

2) SET MAXIMUM TIME IN DATA STATEMEKT MAXSTP; 

EX: DATA MAXSTP /N/ 

( WHERE Ns MAXIMUM NUMBER OF TIME STEPS) 

3) SET INTERVALS FOR PLOTTING BY ENTRIES IN DATA STATEMENT NPLOT; 

EX: DATA NPLOT/ 1,2,8,12,20,40,50,75/. 

( THIS PLOTS OUTPLfT AT THE 1ST, 2ND, 8TH, 1 2TH, 40TH, 50TH, AND . 
75TH PROGRAM TIME STEP) 

4) CHOOSE WHAT OUTPUT DEVICE TO USE, EITHER TEK 618 OR OCMPRES 
(METAFILE OUTPUT) BY GOING TO GRAPHIC DEVICE SECTION OF PROGRAM 
AND COMMENTING OUT WHICH OF TWO DEVICES YOU DONT WANT TO USE. 

5) COMPILE PROGRAM UNDER VS FORTRAN USING THIS COMMAND; 

FORTVS V0R1SPL (OPT (3) 

( NOTE: THIS OPTIMIZES OOOE AT LEVEL 3 ) 

6) GO TO APPROPRIATE TERMINAL AND RUN UNDER DISSPLA. 

A. WHEN THE DISSPLA EXEC PROMPTS YOU ASK FOR 5 CYLINDERS OF 

TEMPORARY DISK STORAGE. 

B. WHEN THE EXEC PROMPTS YOU YOU MAY WISH TO ISSUE: 

FILEDEF 06 DISK V0R10UT LISTING (PERM 

IN ORDER TO REROUTE THE PRINTED OUTPUrT TO A DISK FILE 
NAMED: VCR10UT LISTING A1 . 

0. PROGRAM WILL NOW RUN SENDING OUTPUT TO YOUR DESIGNATED 
DEVICES. 



ADDITIONAL PROGRAM NOTES: 
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c 1 ) 
c 
c 
c 

C 2) 
C 
C 
C 

C 3) 

c 

c 

c 

c 

c 

c 

c 

c 4) 

C 

C 

c 

c 

c 

c 

c 

C 5) 
C 

c 

C 6) 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c««« 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



THIS PROGRAM CCMPOTES DENSITY AM5 VORTICITY IN A (M-2 X N-2) GRID 
BY SOLVING THE BOLTfANTLY SCALED FORM OF THE NAV I ER-STOKES 
EOJATIONS VIA AN "UPWIND DIFFERENCING" FINITE DIFERENCE METHOD. 

THE VELOCITY FIELD- IS COMPLfTED USING A FORM OF THE 0IOT-SAVAT 
LAW AT EACH NDOAL POINT INCLUDING THE CCNTR I BLfT I ON OF THE 
•IMAGE VCRTICIES " ADDED TO INCLUDE THE FREE SURFACE EFFECT. 

THIS VERSION CF THE PROGRAM INCLUDES A OISSPLA RO/TIhC 
TO PLOT THE DENSITY PERTUBATIONS .CONSTANT DENSITY, 

VORTICITY, AND VELOCITY FIELD AT EACH SELECTED TIME STEP. 

THESE SELECTIONS FOR PLOTTING ARE INDICATED BY ENTRIES 
IN THE ARRAY NSTEP. MAXIMUM NUMBER OF TIME STEPS IS INDICATED 
BY AN ENTRY IN THE ARRAY MAXSTEP. 



THE CCMPUTTATIONAL GRID (M-2)X(N-2) IS USED FOR ALL CCMPUTATICNS 
EXCEPT FOR VELOCITY CALLS TO THE BIOT-SAVAT SUBROUTINES. 

TWO ADDITIONAL GRIDS, ONE FOR U AND FOR V, ARE SHIFTED H/2 
UNITS IN THE X AND Y DIRECTION RESPECT I VELY . THE SUBROUTINES 
COMPUTE U AND V VELOCITIES ON THESE GRIDS USING VORTICITY FROM 
OCMPUTATIONAL GRID. THE VELOCITIES ARE THEN AVERAGED ACROSS THE 
COMPUTATIONAL GRID POINTS TO BE USED FOR THE FINITE DIFFERENCING. 

THIS VERSION OF THE PROGRAM IS WRITTEN IN SINGLE PRECISION 
AND IN AN OPTIMIZED FORM CF T>€ ORIGINAL CODE. 

THIS PROGRAM HAS THE OPTION OF SAVING THE RHO AND ZETA ARRAYS 
AND OTHER PROGRAM CONSTANTS IN ORDER TO RESTART THIS SIMULATION 
AT THE TIME A PREVIOUS RUN FINISHED AT. PROGRAM RSVCR1 IS 
SPECIFICALLY WRITTEN TO INPUT THE DATAFILE WRITTEN BY THIS 
PROGRAM AND RESTART THE SIMULATION WITH THE SAME PARAMETERS. 

TO ACCESS THIS SAVE FEATURE SET VARIABLE ISAVE =1 . A FILEDEF 
OF FORM ... FILEDEF 08 DISK <FN> <FT> <FM> ... MUST BE ISSUED 
AT THE APPROPRIATE TIME IN EXECUTION. 

(NOTE: TO DISABLE THIS SAVE FEATURE SET ISAVE =0 ) 



MAJOR VARIABLES USED IN THIS PROGRAM; 



REAL VARIABLES: 

T = NONDIMENSIONAL 
DT = TIME STEP 
H s NONDIMENSIONAL 
RO = NONDIMENSIONAL 
XWIDTH: NONDIMENSIONAL 
ZTMAX = MAXIMUM VORTEX 
SP = NONDIMENSIONAL 
FV = NONDIMENSIONAL 
BVND = NONDIMENSIONAL 



TIME 

GRID LENGTH PARAMETERS DX s DY 
VORTEX CORE RADIUS 
WIDTH CF COMPUTATIONAL AREA 
STRENGTH 

STRATIFICATION PARAMETER 
FROUDE NUMBER 
BRUNT-VALSA FREOJEhO' 



INTEGER VARIABLES; 

M = NUMBER OF NODES IN X DIRECTION 
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CM = NUMBER OF NOOES IN X DIRECTION 

C N s NUMBER OF NOOES IN Y DIRECTION 

C MP s X NODAL POINT OFVORTE>: START POSITION 

C NP s Y NODAL POINT OF VORTEX START POSITION 

0 MM 5 X ARRAY DIMENSION 

c nn' = y array dimension 

0 NFIG s PLOT NUMBER 

0 NSTEP = ITERATION NUMBER 

0 

0 ARRAY VARIABLES; 

0 RHO s DENSITY PERTTJBAT I ONS AT EACH NODE DUE TO DYNAMICS 
0 RHOI s INITAL VALUES OF DENSITY AT EACH NODE DUE TO GRADIENT 

0 RHOT = SUM OF BOTH DENSITY EFFECTS AT EACH NODE 

0 RHICW r NEW VALUE OF RHO AT h€XT TIME STEP 
0 2ETA = VCRTIOITY AT EACH NODE 
C ZTTCW = h€W VALUE OF VORTICITY AT ^€XT TIME STEP 
C U = X OOMPO^eNT OF VELOCITY AT EACH NODAL POINT 

C V s Y OCMPCNENT OF VELOCITY AT EACH NODAL POINT 

C Ua r U VELOCITYS COMPUTED AT BIOT-SAVAT NODAL POINTS 

C V2 = V VELOOITYS COMPUTED AT BIOT-SAVAT NODAL POINTS 

C X s X COORDINATES AT EACH NODAL POINT 

C Y s Y COORDINATES AT EACH NODAL POINT 

C XU s X COORDINATE AT EACH U B. S. NODAL POINT 

C YU = Y COORDINATE AT EACH U B. S. NODAL POINT 

C XV = X COORDINATE AT EACH V B. S. NODAL POINT 

C YV s Y OOCRDINATE AT EACH V B. S. NODAL POINT 

C 

c 

c 

C SUBROUTINES USED IN THIS PROGRAM: 

C 

C SUBROUTIhC OOUT - 

C 

C SUBROUTINE OUTPT - 

C 

C SUBROUTirC CPLOT - 

c 

0 

0 SUBROUTINE ZPLOT - 

C 
C 

C SUBROUTINE DRPLOT- 

C 
C 
C 

C SUBROUTINE ZPL 

C 

C SUBROUTINE VPLOT - CREATES COMPLEX VELOCITY VECTORS AT EACH 

C NODAL POINT AND THEN CALLS THE FOLLOWING 

C SUBROUTINES TO PLOT THEM ON DISSPLA; 

C A) V2PLOT 

0 B) AMIN 

0 0) AMAX 

C 

C SUBROUTINE CON - CALLED BY ZPL TO CONTOUR VORTICITY ON PRINTER. 



OUTPUTS NODE COORDINATES 

OUTPUTS VARIOUS PROGRAM ARRAYS WHEN CALLED 

CREATES ARRAY OF DENSITY AT EACH POINT 

AND CALLS DISSPLA CONTOURING ROUTINE TO PLOT. 

CREATES ARRAY OF VORTICITY AT EACH POINT AND 
CALLS DISSPLA OOUTCURING ROUTINE TO PLOT. 

CREATES ARRAY OF DENSITY PERTUBATlOe AT EACH 
POINT AND CALLS DISSPLA CONTOURING ROUTINE 
TO PLOT. 

OUTPUTS VORTICITY CONTOURS AS A PRINTER PLOT. 
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FUNCTIONS USED IN THIS Pf^OGRAM: 

^FUNCTION VBS - COMPUTES THE V COMPO^€Nr OF VELOCITY USING 
THE BIOT-SAVART LAW. 

FUNCTION UeS - COMPUTES THE U OOMPCtCNT OF VELOCITY USING 
THE BIOT- SAVART LAW. 

FUNCriCN Pt^ - TAKES THE FINITE DIFFERENCE TO STEP THE 
EQUATIONS IN TIME. 



DIMENSION RHOI (25,50) ,U(25,50) ,V(25,50) ,ZETA(25,50) ,RH0(25,50) 
» ,ZTNEW(25.50) ,RHNEW{25, 50) ,NPL0T(8) ,X(25) ,Y(50) 

« ,XU(25) ,YU(50) ,XV(25) ,YV(50) ,U2(25, 50) , V2(25, 50) 

CHARACTER»8 LABEL 
COMMON WK(15000) 

COMMON /PN/ DT, I , J,UF,UFA,UB,U8A,VF,VFA,VB,VBA,VPT,FV 
COMMON /PARM/ H,M,N,PI 
DATA MAXSTP/140/ 

DATA NPLOT/1 ,20,40,60,80, 100, 120, 140/ 

( CALL GRAPHIC OUTPUT DEVICE ) 

CALL CCMPRS 
CALL TEK618 



MM= 25 
NN= 50 

INPUT PROGRAM PARAMETERS 



DT= 2.000E0 
T=0.0 

XWIDTH = 5.0 
R0= 0.09E0 
SP= 0.75 
FV=0.0182 
ZTMIN= 1.0E-8 
M= 24 
N= 44 
I SAVE: 0 

DESIGNATE NOOAL POINTS TO PLACE INITAL VORTEX AT (MP,NP) 
MP=4 
NP=34 

-- COMPUTE ADDITIONAL PARAMETERS 

DENSITY GRADIENT IN Y DIRECTION INITALLY s ORHO 
DRHOI :-SP«SP»FV»FV 

C NONOIMENSIONAL B. V. FREOJENCY r BVND 

BVND= SP«SP»FV»FV 

C MAXIMUM VORTICITY AT ANY POINT = ZTMAX 

ZTMAX=-(FV/(R0»R0)) 
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C GECMETRIC OCNSTANT PI 

PU 4.0E0»ATAN(1 .OeO) 

C GRID LENGTH PARAMETER = H 

H= XWID'm/(M-4) 

CCMPLfTE Wise M AND N PARAMETERS 
MMis M-1 
MPi= M+1 
MM2= M-2 
NPi=N-»-l 
M1i=N-l 
M12s N-2 



INITALIZE CX3DRDINATE SYSTEM ORIGINS 

X(1)= -H 
Y(1)= H 
XU(1)= -1.5E0«H 
TUCDs H 
XV(1)= -H 
YV(i)= 1.5E0»H 

OOMPLfTE COORDINATES 

( NCfTE: COMPUTATIONAL GRID IS IN 4TH GUADRANT ♦X -Y ) 

DO 10 Is 2.MP1 
OO 10 Js 2.NP1 
XU(l)s XLI(I-I) 

YU(J)s TUCJ-D-H 
XV(l)s XV(I-1)+H 
YV(J)s YV(J-1)-H 
X(l)s X(l-1)+H 
10 Y(J)s Y(J-1)-H 

YHT s ASS (Y (N-2)) 

... OLfTPUT INITAL PARAMETERS OF THIS RUN OF PROGRAM 

S*RITE(6. 1275) 

V4?ITE(6, 1500) 

VRITE(6, 1535) OT.H 
\*RITE(6, 1540) XWIDTH.YHT 
V^ITE(6, 1550) SP.FV 
V»RITE(6, 1560) M,N 
WRITE (6, 1570) BVND 

OUTPUT COORDINATES OF NODES 

WRITE(6, 1560) X(MP),Y(NP) 

WRITE (6, 1590) 

CALL COUT(X,Y,MM,NN) 

SPECIFY INITAL V0RTICI7Y AND DENSITY FIELDS 

( FOR GRID SYSTEM X,Y { OCMPUTAT I ONAL GRID) IN THE 4TH QUADRANT ) 
( TERMS ARG1 - AR64 ARE IN QUADRANTS 1 - 4 ) 

ARGMN= 30.0 
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XCENTs X(MP) 

YC£KT= Y(NP) 

DEN= 2.0£0»R0«R0 
2ETA(MP,NP)= ZTMAX 
SO 30 J=1 ,N 
DO 30 1=1 ,M 
TRMUO.OEO 
TRM2=0.0e0 
TT?M3=0.0e0 
TRM4=0.0e0 

APG1= ((X(I)-XCCNT) ii2 ♦ (Y(J) >YC£NT) ««2) /DEN 
IF (ARGl .GT. ARGMN) GO TO 12 
TRMIs ED<P(-AiRG1) 

12 ARG2= ((X(l)+XCENr) i«2 ♦ (Y(J) +YC£Kr) «»2) /DEN 
IF(ARG2 .GT. ARGMN) GO TO 14 
TRM2= EXP(-ARG2) 

14 ARG3= ( (X ( I ) +XCENr) » «2 (Y(J)-YC£NT) *ii2) /DEN 
IF(ARG3 .GT. ARGMN) GO TO 16 
TRM3= EXP(-ARG3) 

16 ARG4= ((X(I)-XCENT) »»2 > (Y(J) -YCENT) »«2) /DEN 
IF(ARG4 .GT. ARGMN) GO TO 1t3 
TRM4= EXP(-ARG4) 

18 IF((J.EQ.NP) .AND. .(I .EQ.MP)) GO TO 20 

2ETA( I , J) =ZTMAX»(-TRM1 +TRM2 -TRM3 *TRM4) 

IF(AeS(2ETA(l , J)) .LT. ZTMIN) ZETA( I , J) =0.0E0 
20 RHOI(l,J)= DRHOI»Y(J) 

RHO(l ,J)= O.OEO 
30 CONTINUE 

VtRITE(6, 1250) 

CALL ajTPT(ZETA, 'ZTA( I , J) ' .MM.hW.O) 

CALL 0LrTPT(RH0l,'R0l(l,J)',MM,hN,0) 

INITALIZE TIME STEP OCUNTZRS 

NFIG=0 

NSTEP=0 



ENTRY POINT OF TIME LOOP^ 

60 CONTINUE 



COMPUTE U VELOCITY FOR ALL NODES IN U VELOCITY GRID 



DO 70 Is 1,MP1 
DO 70 J= 1,N 

70 U2(I,J)= U8S(XU(I) ,YU(J) .MM.NN.X.Y.ZETA) 

COMPUTE V VELOCITY FOR ALL NODES IN V VELOCITY GRID 



DO 100 1= 1,M 
DO 100 J= 1,NP1 

100 V2(I,J)= VBS(XV(I) ,YV(J) ,MM,M^,X,Y,ZETA) 

INTERPOLATE VELOCITIES FROM VELOCITY GRIDS TO NODAL 

POINTS ON COMPUTATIONAL GRID. 



VM I N= 1 . OE-08 
DO 120 I = 1 , M 
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C» 120 J= l.N 

U(I,J)= (U2(I,J) ♦ U2(U1,J))/2.0GO 
V(I,J)= (V2(I,J)+ V2(l ,J*1))/2.0e0 
IF (ABS(U(I,J)) .LT. VMIN) U(I,J)= O.OEO 
_ IF (ABS(V(I,J)) .LT. VMIN) V(l,J)s0.C€0 
120 CCNTINUe 

IF (NSTtP .EQ. 0) CALL OLrTPT(V, W(l ,J) '.MM.M^.l) 
IF (NSTEP .EQ. 0) CALL OLnPT(U, *U( I , J) *,MM,M^,1) 



INCREASE TIME COUNTER BY ONE. 

CHECK ITERATION OOUNTER AND TERMINATE IF NSTEP IS MORE THAN MAX 
COMPUTE VORTICITY ANO DENSITY AT ALL INTERIOR POINTS 
( STORED AS ELEMENTS OF ZT>€W ANO RHh€W ) 



NSTEP=NSTEP*1 

IF (NSTEP .GT. MAXSTP) GO TO 1700 
T= T+DT 
DO 600 lr2.MMl 

DO 600 Js2,NM1 
VPTs V(l , J) 

UF= (U(U1,J)+U(I,J))/2.0E0 
UFA=AfiS(UF) 

Ue= (U(l-1 , J)-KJ(I ,J))/2.0E0 
UBAsABS(UB) 

VFs (V(l , J-1) W(l ,J))/2.0E0 
VFAsABS (VF) 

VBr (V(l ,Jt1)+V(l ,J))/2.0£0 
VBA=ABS(VB) 

2msEW( I , J)= PNEW(ZETA,RHO, MM, M^,-1.0E0, O.OEO, O.OEO) 
RH^€W( I , J)= PNEW(RHO,RHO,MM,NN,O.OEO,O.OEO,BVND) 



— APPLY BOUNDRY CONDITIONS TO RHO ANO ZETA 

. ASSIGN NEW VALUES OF RHO AND ZETA TO RHO ANO ZETA ARRAYS . . . 

DO 700 1= 2.MM1 
DO 700 J= 2, Nil 
ZETAd ,J)= ZTTtW(l ,J) 

RHO(l ,J) :RHr^(l,J) 

CALCULATE RHO AND ZETA VALLES AT I =M-2 ANO J=N-2 BOUNORIES ... 

DO 800 1= 2, MM2 
DO 800 Jr 2.N12 

ZETA(MM2,J)r (ZETA(MM1 ,J)+ ZETA(M-3, J) )/2.0E0 
ZETA(I.NM2)= (ZETAd ,NM1)+ ZETA(I ,N-3) )/2.0E0 
RHO(MM2,J)= (RHO(MMl.J)* RH0(M-3, J) ) /2.0E0 
RH0(l,NM2)r (RHOd.NMl)* RHO( I ,N-3) )/2.0£0 

IF(A8S(ZETA(I , J)) .LT. ZTMIN) ZETA( I , J) rO.OEO 
IF(ABS(RHO(l ,J)) .LT. ZTMIN) RHO( I , J) r 0. OEO 
800 CONTINUE 
C 
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PLOT FLOW PATTERN WHEN NSTEP = VALLE SPECIFIED IN NPLOT 



C 
C 
C 

C THEN RESTART ALC3DRITHM AT LINE 6 FOR A lEV TIME STEP 
• • • • * » ••••••••••••«•••••••••••••••••••••••••••••••••••«••••••••«••••• 

DO 1100 K= 1,0 

IF (NSTEP-NPLOT(K)) 1100, 1200, 1100 

1100 OONTTINUE 
GO TO 60 
1200 hFIG^I^IG+l 
TM= T«FV 

IF (NSTEP .EQ, MAXSTP ) CALL OLfTPr (ZETA , ' ZTA ( I , J) ' , MM , NN , NSTEP) 
IF (NSTEP .EQ. MAXSTP) CALL OLTTFT (RHO, 'RHO( I , J) ' ,MM,NN,NSTEP) 
WRITE(6, 1300) NSTEP, T,TM 
CALL CPLOT(RHO,RHOI ,MM,hW,NSTEP) 

CALL VPLOT(U,V,X,Y,MM,hW) 

CALL DRPLOT(RHO,MM,MM,NSTEP) 

CALL ZPL(2ETA,MM,M'I,XWIDTH,YHT, ITER) 

OLL 2PLOT(2ETA,MM,M^) 

1250 FORMAT (// 10X,'INITAL 2ETA AM) RHO VALUES : *) 

1275 FCRMAT(M') 

1300 FORMAT (/I OX, ' ITERATION NUM: " , 16, 4X, 'Tl ME=" ,F9. 3, 4X, 'T« =",F10.6) 
1400 WRITE (6, 1500) 

1500 FORMAT ( ///) 

1530 FORMAT(' ',5X,'INITAL PROGRAM PARAMETERS:' ) 

1535 F0RMAT('0' ,5X, 'TIME STEP: ' ,F8.4,5X, 'GRID SPACING H = ',F8.4) 
1540 FORMAT (' 0' ,5X, 'NON DIMENSIONAL CELL DIMENSIONS : XDIRECTION s', 

» F8.4,5X,'YDIRECTION s',F0.4) 

1550 F0RMAT('0' ,5X, 'STRATIFICATION PARAMETER s',F0.4,5X, 

» 'FROUDE NUMBER s',F8.4) 

1560 FORMAT (' 0' ,5X, 'NUMBER OF NODES IN X DIRECTION s',|4,5X, 

« ' IN Y DIRECTION s' ,14) 

1570 FCRMAT('0' ,5X, 'NONOIMENSICNAL B. V. FRECUENCY s',F15.9) 

1580 FORMAT (' 0' ,5X, 'X OOOROINATE OF VORTEX START POSITION s',F8.4,5X, 
» 'Y COORDINATE :',F8.4) 

1590 FORMAT (' 0' ,5X, 'OUTPUT OOORDINATES OF ALL OTHER NODES:') 

C 

0 LOOP BAOK TO STARTING POINT OF ITERATION 

GO TO 60 
1700 CALL DONEPL 
C1700 CONTINUE 
C 

C SAVE ARRAYS AND TIME DATA TO RESTART PROGRAM LATER — 

C 

IF (ISAVE .EQ. 0 ) GO TO 1900 
WRITE(0,i») SP,FV,DT,T 
WRITE(0,») M,N,MP,NP 
WRITE(8,») BVND,H,ZTMAX 
DO 1800 1= 1,M 
DO 1800 J: 1,N 

1800 WRITE(8,») 2ETA(I,J), RHO(l,J) 

1900 STOP 
END 

0 



132 



ooo o o oooooooooo ooooo 



SUBROUriNES: 



c 

c 

C 

C 

C 



REAL FUNCTION PNEW(P,Q,MM,hM, A,B,C) 

THIS FUNCTION INTEGRATES THIS EGUATION WITH RESPECT TO TIME 

DP/OT= (-D(UP)/CK -0(VP)/DY ♦ AiOQ/OX + 8»C€L(P) -»CiiV ) $$$ 

UTILIZING AN UPWIND DIFFERENCING SCHEME 

DIMENSION P(MM,lsN) ,Q(MM,M^) 

COMMON /PN/ OT, l,J,UF,UFA,LB,UBA,VF,VFA,VB,VBA,VPT,FV 
COMMON /PARM/ H,M,N, PI 

Pis (UF-UFA) »P(I ,J) ♦(UF+UFA-UB+UBA) iPd , J) -(UB+UBA) »P( I -1 , J) 
P2s (VF-VFA) »P(I , J-1) ♦(VF+VFA-VB+VBA) iPCI , J) -(V8+VBA) »P( I , J+1 ) 
P3s Q(l+1 , J) -0(1-1 ,J) 

P4s P(I+1,J) ♦Pd-l ,J)+P(I ,J-1)+P(I ,J+l)-4.0eO»P(l ,J) 

P5s 2.0E0«H»VPT 

PNEWrPd , J) + (0T/(2.C€0»H)) »(-Pl-P2 +A»P3> (2.0E0/H) iB»P44C»P5/ 

RETURN 

END 



REAL FUNCTION UBS (XU, YU.MM.N^.X, Y.ZETA) 

THIS SUBROUrihE FINDS BOUNORY VALUES OF VELOCITY U 
UTILIZING THE BIOT-SAVAT LAW 

(NOTE: TERMS IN EQUATION ARE EVALUATED CONSECUTIVELY IN EACH 
GUADRANT FRCM CYC TO FOUR ) 



DIMENSION X(MM) ,Y(hN) ,2ETA(MM,hN) 
COMMON /PARM/ H,M,N,PI 

UMINr l.OE-00 
Uls O.OEO 
U2s O.OEO 
U3s O.OEO 
U4s O.OEO 



SIGNs-l.OEO 
DAs H»H 

COMPUTE U CCMPaCNT OF VELOCITY USING THE B-S LAW 

— GUADRANT CNE ~ 

DO 5 Is 1,M 
DO 5 Js 1 ,N 

RSQls((XU-X(D) m.2 4- (YU+Y (J) ) » 12) 

5 Uls UU ((-Y(J)-YU)/(2.0E0»PI»RSQ1))»SIGN«ZETA(l,J)i0A 

C — GUADRANT TWO — 

SIGNS l.OEO 
DO 10 Is i,M 
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DO 10 Js 1 ,N 

RSQ2= ((XU+X(I)) «na ♦(YU+Y(J))»«2) 

10 U2» U2+ ((-Y(J)-YU)/(2.0e0»PI »RSQ2 ))bSIGN»ZETA(I,J)»0A 

C — ajADRANT THREE — 

_§IGN=-1.0E0 
DO 15 1= 1,M 

DO 15 Jr 1,N 

RSQ3r ((XU+X(D) »»2 ♦('nj-Y(J))»»2) 

15 U3r U3+((Y(J)-\U)/(2.0eO»PI»RSQ3))»SIGN»ZETA(l,J)ii0A 

C — OUADRANT FOUR — 

SIGNr 1.0E0 
DO 20 1= 1,M 
DO 20 Jr 1,N 

RSQ4r ((XU-X(D) »»2 ♦ (YU-Y(J) ) *«2) 

20 U4 r U4 + ((Y(J)-YU)/(2.0EO*PI »RS04)) »SIGN»ZETA(I , J) »DA 
LfTEMPr U1HJ2 HJ3+ U4 
IF (A8S(UTEMP) .LT.UMIN) UTEMPrO.OEO 
UBSr IJTEMP 
RETIJRN 
END 
S' 



REAL FUNCTION VBS(XV,YV,MM,I^,X, Y.ZETA) 

THIS SUBROLfTIhEI FINDS BCUNDRY VALUES OF VELOCITY V 
UTILIZING THE BIOT-SAVAT LAW 

(NOTE : TERMS IN THE EOUATICN ARE EVALUATED OONSEOLfT I VELY IN EACH 
OUADRANT FROM ONE TO FOUR) 



DIMENSION X(MM) ,Y(NN) .2ETA(MM,hM) 
OOMMON/PARM/ H,M,N,PI 



VMINr 1.0E-08 
VI = O.OEO 
V2r O.OEO 
V3r O.OEO 
V4r O.OEO 



SIGNr-1 .CEO 
DAr H«H 

COMPUTE V COMPONENT OF VELOCITY USING THE B-S LAW 

~ OUADRANT ONE — 

DO 5 Ir 1 ,M 
DO 5 Jr 1,N 

RSQ1r((XV-X(l)) »ii2 ♦ (YV+Y(J) ) »»2) 

5 V1= V1+ ((XV-X(I))/(2.E0«PI*RSQ1))»SIGN«2ETA(I,J)»DA 

0 — OUADRANT TWO ~ 

S I GN= 1 . OEO 
DO 10 Ir 1,M 
DO 10 Jr 1 ,N 

RSQ2= ((XV+X(I))»»2 +(YV>Y(J) ) «»2) 
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10 WZz ((XV+X(l))/(2.0e0»PliiRSQ2))»SIGN»2ETA(l,J)»0A 
C — ajADRAfiT TVIREE ~ 

SIGNr-I.OEO 
DO 15 1= 1,M 
_ DO 15 J= 1,N 

RSQ3= ((XV+X(D) ■■a >(YV-Y(J))«»2) 

15 V3s V3*((XV+X^I))/(2.0e0iiPliiRSQ3))»SI<jN»2ETA(l,J)»DA 
C — GUADRAKT FOUR — 

SIGNS i.oeo 

DO 20 Is 1 ,H 
DO 20 Js 1 ,N 

RSQ4s ((XV-X(D) »»2 ♦ (YV-Y(J)) »«2) 

20 V4 s V4 ♦ ((XV-X(I))/(2.0E0»PI»RSQ4))«SIGN»ZETA(I,J)»DA 
VTEMP = V1+V2 +V3+ V4 
IF(A8S(V7EMP) .LT.VMIN) VTEMPsO.OEO 
VeSs VTEMP 
RETURN 
END 



SUBROUTINE OOUT(X,Y,MM,hM) 



REAL X(MM) ,Y(W) ,H,PI 
OOMMON/PARM/H , M , N , P I 
VvRITE(6,50) 

WRITE(6,«) ' ARRAY X' 

V»RITE(6, ») 

VRITE(6,») (X(l), Isi.M) 
V»RITE(6, ») 

WRITE(6, «) 

V*RITE(6,») ' ARRAY Y 

WRITE (6,») 

WRITE(6, ») (Y(J) , Js1 ,N) 

WRITE(6, ») 

W?ITE(6, «) 

50 FORMAT (///) 

RETURN 

END 



SUBROUTINE OUTPT(A,LBL,MM,NN,KAL) 

TVIIS SUBROUTINE PRINTS OUT AN (N<M) ARRAY WITH A LABEL 
IN A RECTANGULAR ( 2,M-2)X(2,N-2) GRID 

DIMENSION A(MM,M^) 

COMMON/PARM/ H,M,N,PI 
CHARACTER »8 LBL 
NM2= N-2 
MM2s M-2 
WRITE(6,50) 

WRITE(6, ») 
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CLfTPOr or ARRAY: 



,LBL, 



ITER MJM: ' ,KAL 



V«ITE(6,») ' 

WRITE(6, M) 

DO 10 J» 2, Mia 
UR1TE(6, ») 

^ITE(6,») (A(I,J), 1=2, MM2) 
10 OOKTINUE 
50 F0RMAT.(M^) 

RETURN 

END 



SUBROUTIhC CPLCrr(RHO,RHOI ,MM,M4,NSTEP) 

COMMON WK( 15000) 

COMMON /FARM/ H,M,N,PI 
REAL W(41 ,41) 

REAL RHO(MM,M0 ,RHOI (MM.M4) ,H,PI ,RHCrT(25, 50) 

DATA I SCALE/4HSCAL/ 

. . .CREATE ARRAY OF DENSITY OHANGES AT EAOH NODE (DRHO) 

MM2= M-2 
NM2= N-2 
MW= 41 
hW= 41 

DO 10 1= 2, MM2 
DO 10 J= 2,NM2 
JJ= (M12>1)-J 
11= l>19 

RHOTd , J)=RHO(l , J) ♦RHOI (I.J) 

W(l I , JJ)= RHOT(l ,J) 

10 CONTINUE 

DO 20 1= 3, MM2 
11= 23-1 
DO 20 J= 2,NM2 
JJ= (MI2^-1) -J 
20 W(l I , JJ)= RHOTd , J) 

CALL OL/rPT(RHOT, 'RHT(I,J)',MM,NN,NSTEP) 

CALL SUeROLfTINE CCNTR TO PLOT THESE PERTUBATIONS 

CALL RESET (3HALL) 

CALL PAGE (8. 7, 11.2) 

CALL PHYS0R(2.0,3.00) 

CALL AREA2D (5. 0,6.0) 

CALL HEIGHT(.2 ) 

CALL INTAXS 

CALL XNAME ('X/BO',4) 

CALL YNAME ('Y/BO',4) 

CALL XTICKS(2) 

CALL YTICKS(2) 

CALL GRAF (-5.0, 1 . , 5. , 10. , 1 . ,0. ) 

CALL FRAME 

CALL 800M0N (15000) 

CALL 00NMIN(2.0) 



I 
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CALL OCMANG (30.) 

CALL OWAK (W.MW.NW, I SCALE) 

CALL HEIGHT(.08) 

CALL OONLIN (0, 5HS0L 1 0, 'LABELS' , 2, 5) 

^L COM. IN (1 ,4HDASH,8HN0LA8£LS, 1 ,3) 

CALL RASPLN(0.25) 

CALL OCMTUR ( 2 , 6KjABELS , 4HDRAW) 

CALL HEIGHT(.2) 

CALL OOHPLX 

CALL RESET ('HEIGHT') 

CALL RESET ('COMPLX') 

CALL ENOPL (0) 

RETURN 

END 

C 

C 

C 

SUBROLTfll^ ZPLOT(ZETA,MM,W) 

REAL 2ETA(MM,I^) ,H,PI 
OCMMON WK(1500) 

COIrCN /PARM/ H,M,N,P1 
DIMENSION W(41,41) 

DATA I SCALE/4HSCAL/ 

C 

C. ..CREATE ARRAY FOR PLOTTING OF ZETA IN ACTUAL (XMPUTAT I ONAL AREA.. 
0 

MWs 41 
N*/= 41 
MM2= M-2 
NM2=N-2 

DO 10 Ir 2, MM2 
lls U19 
DO 10 Jz 2,1^2 
JJ= (NM2+1)-J 
to W(l I , JJ)=2ETA(I ,J) 

DO 20 1 1 3, MM2 
I U 23-1 
DO 20 J= 2.NM2 
JJz (NM2>1) - J 

20 W(l I , JJ) = -ZETAd ,J) 

0 

0 CALL SUBROLfTIhC OCNTR TO PLOT VORTICITY 

0 

0 

CALL RESET (3HALL) 

CALL PAGE(8.7, 11.2) 

CALL PHYSCft (2. 0,3.00) 

CALL AREA2D (5. 0,6.0) 

CALL HEIGHT(.2 ) 

CALL INTAXS 

CALL XNAME ('X/BO',4) 

CALL YNAME ('Y/BO',4) 

CALL XTICKS(2) 

CALL YTICKS(2) 

CALL GRAF (-5.0, 1 . , 5. , 10. , 1 . ,0. ) 
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CALL FRAME 

CALL BOCMON (15000) 

CALL OONANG (60.) 

CALL 0t>#ilN(2.0) 

_CALL OOMIAK (W, MW, M*/, I SCALE) 

CALL HEIGKT(.Oe) 

CALL COLIN (0,5HSOLID,'LABELS',2,5) 

CALL COLIN (1,4H0ASH,* LABELS',!, 3) 

CALL RASPLN(0.25) 

CALL COrrUR (2,6HLA8ELS,4HDRAW) 

CALL HEIGKT(.2) 

CALL CCMPLX 

CALL RESET ('HEIGHT') 

CALL RESET ('CCMPLX') 

CALL ENOPL (0) 

RETURN 

END 

C 

C 

C" 

C 

SUBRCLfTINE DRPLOT(RHO, MM, NN, NSTEP) 

COMMON WK( 15000) 

COMMON /FARM/ H,M,N,PI 
REAL W(41,41) 

REAL RHO(MM,NN) ,H,PI 
DATA ISCALE/4HSCAL/ 

C 

C. ..CREATE ARRAY CF DENSITY CHANGES AT EACH NOOE (DRHO) . 
C 

MM2s M-2 
NM2= N-2 
MW= 41 
NW= 41 

DO 10 1= 2, MM2 
DO 10 Jr 2,^t12 
JJr (NM2+D-J 
Nr I-H9 

W(l I , JJ)r RHO(l ,J) 

10 CONTINUE 

DO 20 I r 3, MM2 
Hr 23-1 
DO 20 Jr 2,NM2 
JJr (NM2>1) -J 
20 W(l I , JJ)r RHO(l , J) 

C 

C CALL SUBROUrrirL OONTR TO PLOT THESE PERTUBATIONS. ; 

C 

CALL RESET (3HALL) 

CALL PAG£(8.7, 1 1 .2) 

CALL PHYSOR (2.0, 3.00) 

CALL AREA20 (5.0, 6.0) 

CALL HEIGHT(.2 ) 

CALL I NTAXS 

CALL XNAME ('X/BO' ,4) 
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CALL 'r>Wl£ ('Y/BO',4) 

CALL XTIO:S(2) 

CALL YTICKS(2) 

CALL GRAF (-5.0, 1 . , 5. . 10. , 1 . ,0. ) 

C*J-L FRAME 

CALL BOaCN (15000) 

CALL OMIN(2.0) 

CALL OCNANG (20.) 

CALL OOtWIC (W.MW.NW, ISCALE) 

CALL HEIGHT(.08) 

CALL OCM.IN (0,5HS0LID,/LABELS',2,5) 

CALL cm. IN (1 ,4H0ASH,8WCLA8ELS, 1 ,3) 

CALL RASPLN(0.25) 

CALL OONTUR (2, 6HLABELS.4H0RAW) 

CALL HEIGHT(.2) 

CALL OOMPLX 

CALL RESET ('HEIGHT") 

CALL RESET ('COMPLX') 

CALL ENOPL (0) 

RETURN 

END 

C 

C====:===========rzzzz=z======zzzrrr=======z ======== ====:=====:=:====== 

C 

SUBROUTINE 00N(22,XRANGE,YRANGE,N<.H,MM.M4,ZMAX.2MIN) 

C 

C TWIS SUBROUTIIE GENERATES A COMTCUR PLOT OF ZETA IN A RECTANGULAR 

C . (X>1AIN OF SIZE; YRANGE » XRANGE ON THE PRINTER. 

C 

CHARACTER*! SYMBOL (8) .GRAPH ( 100) 

REAL LX.LY 

DIMENSION Z2(MM,M^) ,VALUE(7) 

DATA SYMBOL / 'X' . 'A' , '♦','1',' '/ 

VALUE (1)= ZMIN 
VALUE(2)= ZMIN*. 75 
VALUE (3) = ZMIN*. 5 
VALUE (4) = 0.0 
VALUE (6) = ZMAX».75 
VALUE (7) = ZMAX 
VALUE(5)=2MAX*.45 
NY= 0. 8 *Nt<* YRANGE ARANIGE 
DELXs XRANG£/N< 

DELYs YRANGE/NY 
V^ITE(6,50) 

V4RITE(6,*) 'N<=',Nt<, 'NTz'. NT, 'XRANIGEr', XRANGE, 

■ 'YRANGE=', YRANGE 

V^ITE(6,») 

DO 11 1= 1.N< 

11 GRAPH(I)= 'T' 

\WRITE(6,6) (GRAPH(I), 1= 1.N«) • 

DO 5 JSYMBL = 1 .NY 
Y= (JSYMBL-0.5) *OELY 
J= 1+Y/H 
LY= Y-(J-1)*H 
HMLY= H-LY 
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DO 4 ISYMSL = t.NOC 
X= (ISYMBL-0.5) »C€LX 
I = 1 ♦X/H 
LX= X-(l-i) nH 
_ HMLXs H- LX 
AU HMLXnHMLY 
A2= HMLXuLY 
A3s LXuLY 
A4= LX«HMLY 

22Cs (A1»Z2(I,J) M2iZ2(l,J+1) M3»Z2( I *1 , J*1 ) 
■ ♦A4«Z2(I*1 ,J))/H»i2 

. . . DETERMIhE THE VALUE OF MRANGE BASED ON 22C 



. ... PRINT CUT CHARACTERS IN ARRAY GRAPH BY ROW 



\ 

2 

3 

4 

5 

6 



50 

50 



DO 2 K= 1,7 

IF (Z2C-VALUE(K)) 1, 1, 2 
M?ANGE = K 
00 TO 3 
CONTINUE 
M?ANGE= 0 

GRAPH ( ISYHBL)= SYH90L (NRANGE) 
CONTINUE 

VRITE(6,6) (GRAPH(I), 1=1, N<) 
FORMATC ',10X,100At ) 

FCRMAT(' 1 ') 

DO 60 1= 1,N< 



GRAPH ( I ) 
WRITE(6,6) 
V^ITE(6, ») 
WRITE(6, ») 
« ,ZMIN 
VRITE(6, ») 
VRITE(6, ») 
VRITE(6, ») 
M?ITE(6 , ») 
« , 2MAX 
VyRITE(6, ») 
RETURN 
END 



'B' 

(GRAPH(I) ,1 = 1, NX) 

'MAXIMUM VALUE OF NEGATINE VORTICITY = SYMBOL (X) 

' 75'/ OF MAXIMUM NEGATIVE VCRTICITY = SYMBOL (A) ' 
' 50'/. OF MAXIMUM NEGATIVE VCRTICITY = SYMBOL (-) ' 

' MAXIMUM VALUE OF POSITIVE VORTICITY = SYMBOL (1) 

' 75’/. OF MAXIMUM POSITIVE VCRTICITY = SYMBOL (+)' 



SUBROUTINE 2PL(ZETA,MM,NN,XWIDTH,YHT, ITER) 

COMMON /PARM/ H,M,N,PI 
DIMENSION 22(25,50) ,2ETA (MM, NN) 

. . .CREATE ARRAY FOR PLOTTING OF ZETA IN ACTUAL COMPUTATIONAL AREA. . 



MM3= M-3 
NM3:N-3 
DO 10 I = 1 ,MM3 
DO 10 J= 1 ,NM3 

10 22(1 ,J)= 2ETA(I f1 ,J+1) 

0. .FIND THE MAXIMUM POSITIVE VORTICITY IN ARRAY 22(MM,NN) 
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c 



ZMINs 0.0 
2MAX= 0.0 
DO aO U 1.MM3 

__ DO 20 J= 1.NM3 

IF(22(I,J) .GT. ZMAX) 2MA>^=Z2 ( I , J) 

20 OOfTINUE 

. FIND THE MAXIMUM hCGATIVE VCRTICITY IN ARRAY Z2(MM,W) 

DO 40 Is 1.MM3 
DO 40 Js 1,M13 

IF( Z2(I,J) .LT. 2MIN) ZMINs Z2(I,J) 

0 CONTINUE 

IF (ZMAX .LT. 1.0E-4 ) ZMAX=1.0 
CALL SUBROLfril^ CON TO PLOT VORTICITY. 

CALL OONCZa.XWIDTH.YHT, 100,H.MM,I^,ZMAX,ZMIN) 

RETURN 

END 



SUBROLfTINE VPLOT(U, V,X, Y,MM,I^) 

COMMON /FARM/ H.M.N.PI 
COMPLEX ZV (21 ,41) , CV(21,41) 

DIMENSION U(MM,hM) , V(MM,M'0, X(MM), Y(M4) 
MM2s M-2 
I^2s N-2 



FORM COMPLEX U, V, AND Z ARRAYS FOR OJADRANTS 3 AND 4 

DO 10 Is 2, MM2 
DO 10 Js 2.M12 
I Is 1-1 

JJs (NM2+1) -J 
ZV(ll,JJ)s CMPLX(X(I) ,Y(J)) 

CV(ll,JJ)s CMPLX(U(I , J) ,V(I ,J)) 



.. PLOT VECTORS USING SUBROUTINE V2PL0T 
CALL V2PLOT(ZV,CV,21 ,41) 

RETURN 

END 



SUBROLTriNE V2PL0T(Z,C,M,N) 



» SUBROUTINE V2PL0T TO PLOT MXN OOMPLX VELOCIT FIELD C(M,N) 
» OVER A CCMPLX DOMAIN Z(M,N) 



COMPLEX Z(21 ,41) ,C(21,41) ,CMAX 
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REAL X(50) ,Y(50) ,YX(50) ,^(50) ,X<(50) ,XN(50) ,X2(50) ,Y2(50) 
INTCGER PMAX.PMIN 
E>rrERNAL AMAX 
EXTERNAL AMIN 
C 

DO 20 J =1,N 
DO 10 I si ,M 

X(l) sABS(REAL(C(l,J))) 

Y(l) sASS(AIMA6(C(l,J))) 

10 OOKTINUe 

CALL AMAX(Y,M,YMAX,PMAX) 

CALL AMAX(X,M,XMAX,PMAX) 

YX(J) s YMAX 
XX (J) s XMAX 

. CALL AMIN(Y,M,YMIN,PMIN) 

CALL AMIN(X,M,XMIN,PMIN) 

YN(J) sYMIN 
XN(J) sXMlN 
20 OC>(TlNUe 

CALL AMAX(YX,N.YMAX,PMAX) 

C CALL AMIN(YN,N,YMIN,PMIN) 

CALL AMAX (XX, N, XMAX, PMAX) 

CALL AMIN(XN,N,XMINjPMIN) 

CMAX s CMPLX (XMAX, YMAX) 

UMAX sSQRT( (XMAX) i»2+ (YMAX) ihi2) 

DO 24 I si ,M 
DO 22 Js1,N 

C(l ,J) s C(I,J)AJMAX 
22 CCNTINUE 

24 COfTINUE 

C 

DO 40 J s 1 , N 
DO 30 I s 1 , M 

X(l) sREAL(Z(l , J)) 

Y(l) sAIMAG(Z(l,J)) 

30 CONTINUE 

CALL AMAX (Y,M, YMAX, PMAX) 

CALL AMAX (X,M, XMAX, PMAX) 

YX(J) s YMAX 
XX (J) s XMAX 
CALL AMIN(Y,M,YMIN,PMIN) 

CALL AMIN(X,M,XMIN,PMIN) 

YN(J) sYMIN 
XN(J) sXMlN 
40 CONTINUE 

CALL AMAX (YX,N, YMAX, PMAX) 

CALL AMIN(YN,N,YMIN,PMIN) 

CALL AMAX (XX, N, XMAX, PMAX) 

CALL AMIN(XN,N,XMIN,PMIN) 

C 

CALL RESET (3HALL) 

CALL PAGE (8. 7, 11,2) 

CALL PHYS0R(2, 0,3.0) 

CALL AREA2D(5. ,6. ) 
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CALL HEIGHT(.2) 

CALL INTAXS 

CALL XNAME('X/BO' ,4) 

CALL YIWC('Y/0O' ,4) 

O^LL XTICKS(2) 

CALL yric*:s(2) 

CALL GRAF(0.0,0.5,5.0,-10.0, 1.0.0.0) #« 

C 

DO 60 J=i ,N 
DO 50 1=1 ,M 

X(l) = REAL(Zd.J)) 

X2(l) = X(l) ♦REALCCd ,J)) 

Y(l)= AIMAG(2d ,J)) 

Y2(l) = Y(l)MIMAG(Cd, J)) 

50 ocr/riNue 

CALL MARKER (3) 

CALL SCLPtC(.OOi) 

CALL CURVE 0<,Y,M,-1) 

DO 54 I =1,M 
XI 1=X(I) 

X12iX2d) 

Y1 1=Y(I) 

Y12=Y2(I) 

IF(X12.EQ.0. .AN0.Y12.EQ.0) GO TO 54 
CALL RLVEC (XII ,Y1 1 ,X12,Y12, 1 121 ) 

54 CONTINUE 

60 CONTINUE 

CALL RESET ('HEIGHT') 

CALL ENDPL(O) 

RETURN 

END 

0 

0 

SUBROUTINE AMAX (Y,N, YMAX, PMAX) 

0 

0 » SUeROUTIh€ TO OCNRJTE THE LARGE ELELEMENT 

0 a IN A GIVEN RCW OR COLUMN IN ARRAY 'A' . 

0 aaaaaaaaiaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

0 

0 ««a VARIABLE DECIjARATION «aa 

0 REAL AMAX 

INTEGER PMAX 
DIMENSION Y(N) 

0 

YMAX=Y(1) 

DO 5100 1=2, N 

IF(Yd) .LE.YMAX) GO TO 5000 
YMAX=Y(I) 

PMAX=I 

5000 CONTINUE 

5100 CONTINUE 

0 MAXRCWrPOS 

RETURN 
END 
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c 

c 

c 



SUBROLfTIKE AMIN(Y,N,YMIN,PMIN) 



C 

C 

C 

C 

C 

C 

C 

C 



■«« VARIABLE DECLARATION 
REAL AMIN 
INTEGER PMIN 
DIMENSION Y(N) 



» SUeROUTII^ TO OCMF=UTE THE SMALL ELEMENT 
« IN A GIVEN ROW OR COLUMN IN ARRAY 'A' . 



0 



YMIN=Y(1) 

DO 6100 1=2, N 

IF(Y(I) .GE.YMIN) GO TO 6000 
YMIN=Y(I) 

PMINrI 

6000 CONTINUE f 

6100 CONTINUE 



RETURN 

END 



EOF; 
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